Climate: Degree Heating Weeks

Introduction

This script contains the methods used to wrangle, analyse and present degree heating week (DHW) data in the Northern Three regions. For a guide on downloading DHW data refer to the README document for the Spatial Analysis GitHub repo. Note that DHW data should be downloaded automatically by this script.

Degree Heating Weeks Data is predominantly used within the climate section of the technical report as a proxy for coral bleaching. The description of a degree heating week is a bit confusing, but as defined by NOAA: “DHWs show the accumulated heat stress experienced by corals in the prior three months, it is a cumulative measure of both intensity and duration of heat stress”.

Temperatures exceeding 1°C above the usual summertime maximum are sufficient to cause stress, including bleaching, and are the basis of a degree heating week. For example, a DHW of 2 is equivalent to one week of Hot Spot values persistently at 2°C, or two weeks of Hot Spot values persistently at 1°C above usual summertime maximum temperatures. Values over 4 have been shown to cause significant coral bleaching, and values over 8 have caused severe bleaching and significant mortality.

The main objectives of this script are to:

• Automatically download all required DHW data
• Bin values into a clear and conscience legend
• Create a map of the current financial years’ degree heating weeks
• Create a map of the past 5 years of degree heating weeks

Short History

Note

October 2022: This data is currently provided free of charge by Angus Thompson at AIMS. He will provide a 5 year history along with other data such as coral. However, it is important that we understand and can effectively mimic his data should issues arise. This also allows us further customization of the end product without bothering Angus.

Note

February 2023: We have been cleared to reproduce the same analysis in house, and are currently looking at changing the styling to match other in house analyses.

Script Set Up

This script requires multiple core spatial packages, each of which is loaded below. Key variables and paths are also created.

#use pacman function to load and install (if required) all other packages
pacman::p_load(tidyverse, glue, here, janitor, sf, tmap, exactextractr, RColorBrewer, stars, ggplot2)

#install/load the custom package
#pak::pak("add-am/RcTools")
library(RcTools)

#define the target year
script_fyear <- 2025

#build the data path and output path variables
data_path <- here("data/dhw/")
output_path <- here(glue("outputs/dhw/{script_fyear}/"))

#and create some additional sub folders
dir.create(glue("{data_path}/raw/"), recursive = TRUE)
dir.create(glue("{data_path}/processed/"))
dir.create(glue("{output_path}/plots/"), recursive = TRUE)
dir.create(glue("{output_path}/maps/"))

Load Data

Now the script is set up we need to load in all of the required datasets. This will be broken into two segments:

• Spatial data specific to the N3 region - such as the region, basin, and sub basin boundaries.
• Degree Heating Week data

Spatial Data

The northern three boundaries are built using a custom function or reloaded from file.

#if the data is on file, load it, otherwise built it
if (file.exists("data/n3_region.gpkg")){
  n3_region <- st_read("data/n3_region.gpkg")
} else {
  n3_region <- build_n3_region()
}

The dataset is then divided into sections for ease of use later on.

#cut down the n3 dataset to just the marine region
n3_marine <- n3_region |> 
  filter(Environment == "Marine",
         BasinOrZone != "Burdekin Marine",
         Region != "Burdekin") |> 
  rename(Zone = BasinOrZone) |> 
  mutate(Region = case_when(Region == "Burdekin" ~ "Dry Tropics",
                           T ~ Region)) |> 
  group_by(Region) |> summarise(geom = st_union(geom)) |> 
  ungroup() |> st_cast()

#cut down to everything but the marine region
n3_basins <- n3_region |>
  filter(Environment != "Marine") |> 
  rename(Basin = BasinOrZone) |> group_by(Region, Basin) |> summarise(geom = st_union(geom)) |> 
  ungroup() |> st_cast()

#read in reef data and update crs
reefs <- get(data("gbr_feat", package = "gisaimsr")) |> 
  name_cleaning() |> 
  filter(FeatName == "Reef") |> 
  st_transform("EPSG:7855")

NoteWT Request

A special request from the WT team is to use their marine boundary, rather than the marine boundary as defined by the NRM groups. Thus below we need to replace to NRM WT boundary with the custom version. Note that this outline was provided by the WT team and has no online equivalent.

#
#read in the custom WT outline and edit to match the main dataset
wt_custom <- st_read(here("data/historical_wt_boundary.gpkg")) |> 
  name_cleaning() |> 
  select(Name) |> 
  rename(Region = Name)

#remove old WT and add new
n3_marine <- n3_marine |> 
  filter(Region != "Wet Tropics") |> 
  bind_rows(wt_custom)

Degree Heating Week Data

Load in the degree heating week data from NOAA’s website. A custom function has been written that will automatically download the data from the the online servers, then crop the data. If the original or cropped version already exists on your local computer it will not re-download.

#define the place where full files should be saved
full_file_loc <- glue("{data_path}/raw/")

#define the place where cropped files should be saved
cropped_file_loc <- glue("{data_path}/processed/")

#define the area in which you want to crop data to
n3_region <- n3_region

#download, save, crop, and build the data
n3_dhw <- dhw_extract(full_file_loc, cropped_file_loc, n3_region)

#update the crs of the data to match the n3 region geometry
n3_dhw <- st_warp(n3_dhw, crs = "EPSG:7855")

Edit Data

Although we download all DHW data from 1986 to present, currently we are only focused on the most recent 5 years of data. There is no long-term analysis of the data (yet). Below we filter by year to only get the most recent 5 years.

#find the index for the layer that matches the target year
target_layer <- which(paste0(script_fyear, "-06-01") == st_get_dimension_values(n3_dhw, "time"))

#extract that layer plus the four previous ones
n3_dhw_5y <- n3_dhw[,,,(target_layer-4):target_layer]

#use the high res crop function to crop and increase resolution at the same time
n3_dhw_5y <- nc_high_res_crop(n3_dhw_5y, n3_marine, 5)

DHW data is provided with a scale from 0 to ~25, this is a bit excessive for our requirements and makes for a vague and messy legend. Angus Thompson (AIMS) has developed a binned legend to improve clarity of each group, these bins are as follows:

• “Low risk: 0 - 2 DHW”
• “Warning: 2 - 4 DHW”
• “Possible: 4 - 6 DHW”
• “Probable: 6 - 8 DHW”
• “Highly Likely: >8 DHW”

With the legend titled “Bleaching Risk”. The code below bins values in the raster.

Note

NOAA have since updated the scale used for DHW as well, it is pending if we should adopt the same scale.

#extract the underlying array
vals <- n3_dhw_5y[[1]]

#put data into bins
vals_binned <- array(
  findInterval(vals, c(-Inf, 2, 4, 6, 8, Inf), rightmost.closed = TRUE), 
  dim = dim(vals), 
  dimnames = dimnames(vals))

#put the binned values back into the netcdf
n3_dhw_5y[[1]] <- vals_binned

Visualise Data

The two key visualizations this script creates area current year map and a past 5 years map. However before mapping, the legend labels and colours need to be set up.

#dhw colours as per categories and colour blind friendly brewer
dhw_cols <- c("#2C7BB6", "#ABD9E9","#FFFFBF","#FDAE61","#D7191C")

#dhw labels for each bin
dhw_lab <- c("Low likelihood of bleaching (0 - 2 DHW)",
             "Bleaching warning likely (2 - 4 DHW)",
             "Bleaching possible (4 - 6 DHW)",
             "Bleaching probable (6 - 8 DHW)",
             "Severe bleaching likely (>8 DHW)")

#create a vector of the three n3 regions
n3_marine_names <- unique(n3_marine$Region)

Current Finacial Year

First we can plot a single year of data to get an idea of how things look.

#using unique regions
for (i in n3_marine_names) {
  
  #filter by region
  region_basins <- n3_marine |> filter(Region == i)
  
  #get the associated basins
  basins <- n3_basins |> filter(Region == i)
  
  #crop to the specific region and take a single year of data
  single_year_gbr <- st_crop(n3_dhw_5y[,,,5], region_basins)
  
  #mask reefs to the region as well
  target_reefs <- st_intersection(reefs, region_basins)
  
  #plot
  map <- tm_shape(single_year_gbr) +
    tm_raster(col.scale = tm_scale_intervals(values = dhw_cols,
                                             breaks = c(1:6),
                                             labels = dhw_lab),
              col.legend = tm_legend(reverse = T,
                                     title = "Coral bleaching likelihood \n and number of DHW's")) +
    tm_shape(qld) +
    tm_polygons(fill = "grey80", 
                col = "black") +
    tm_shape(region_basins, is.main = T) +
    tm_borders(col = "black") +
    tm_shape(basins) +
    tm_polygons(fill = "grey90", 
                col = "black") +
    tm_text("Basin", 
            xmod = -1.8, 
            ymod = 0.1, 
            size = 0.7,
            options = opt_tm_text(shadow = T)) +
    tm_shape(reefs) +
    tm_borders(col = "grey60",
               col_alpha = 0.4,
               fill = "grey60",
               fill_alpha = 0.2) +
    tm_shape(target_reefs) +
    tm_borders(col = "black", 
               fill = "grey60",
               fill_alpha = 0.2) +
    tm_layout(legend.frame = T, 
              legend.bg.color = "White", 
              legend.text.size = 0.7, 
              legend.position = c("right", "bottom"),
              asp = 1.1)
  
  #edit variable name for better save path
  i <- tolower(gsub(" ", "-", i))
  
  #save to unique variable for later
  assign(glue("mean_map_{i}"), map)
  
  #save the map as a png
  tmap_save(map, filename = glue("{output_path}/maps/{i}_dhw.png"))
  
}

See below for an example of a one year map.

map

Most Recent 5 Years

Then we can create the 5 year map, please note that this was one of the previous outputs given by the AIMS team, which is the main reason this extra map is created. To make this map we can put each year side by side on one map to create a comparison of the currently targeted year and the 4 years proceeding it.

#using unique regions
for (i in n3_marine_names) {
  
  #filter all basins by region
  region_basins <- n3_marine |> filter(Region == i)
  
  #get the associated basins
  basins <- n3_basins |> filter(Region == i)

  #mask to just the region (keep all years)
  all_year_region <- st_crop(n3_dhw_5y, region_basins)
    
  #plot
  facet_map <- tm_shape(all_year_region) +
    tm_raster(col.scale = tm_scale_intervals(values = dhw_cols, 
                                             breaks = c(1:6),
                                             labels = dhw_lab),
              col.free = FALSE,
              col.legend = tm_legend(title = "Coral bleaching likelihood and number of DHW's")) +
    tm_shape(qld) +
    tm_polygons(fill = "grey80",
                col = "black") +
    tm_shape(region_basins, is.main = TRUE) +
    tm_borders(col = "black") +
    tm_shape(basins) +
    tm_polygons(fill = "grey90", 
                col = "black") +
    tm_shape(reefs) +
    tm_borders(fill = "grey60", 
               fill_alpha = 0.6) +
    tm_layout(panel.labels = year(st_get_dimension_values(n3_dhw_5y, "time"))) +
    tm_facets_hstack()

  #edit variable name for better save path
  i_lower <- tolower(gsub(" ", "-", i))
 
  #save the map as a png
  tmap_save(facet_map, filename = glue("{output_path}/maps/{i_lower}_dhw_fyear-{script_fyear}-to-{script_fyear-4}.png"))
  
}

Here is how that facet map looks.

facet_map

Save Data

For the purposes of our report we will also need to save a summary table with some key statistics.

#create a function that extract the min, 25%, mean, 75%, and max, and adds the target region and year
process_region_year <- function(region_name, year_index) {
  
  #what region are we looking at?
  region_basins <- n3_marine |> 
    filter(Region == region_name)
  
  #mask to the specific region and year
  single_year_region <- st_crop(n3_dhw_5y[,,,year_index], region_basins)
  
  #extract the values for the region and calculate quantiles and create a data frame
  quantiles <- quantile(single_year_region[[1]], na.rm = TRUE)
  
  #convert into a dataframe, this is also the output
  data.frame(Region = region_name,
             Year = str_split(st_get_dimension_values(n3_dhw_5y, "time")[year_index], "-")[[1]][1],
             Min = quantiles[1],
             Per25 = quantiles[2],
             Mean = quantiles[3],
             Per75 = quantiles[4],
             Max = quantiles[5],
             row.names = NULL)
}

#create a data frame that has pairs for each region and year combination
region_year_pairs <- expand.grid(Region = n3_marine_names, Year = 1:5)

#run the custom function over every pair in the region year dataset, then bind the outputs together
summary_statistics <- map2(region_year_pairs$Region, region_year_pairs$Year, process_region_year) |> bind_rows()

#save the table
write_csv(summary_statistics, glue("{output_path}/dhw_summary.csv"))